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ABSTRACT 

Although the Earth's orbit is never far from circular, terrestrial planets around other stars 
might experience substantial changes in eccentricity. Eccentricity variations could lead to climate 
changes, including possible "phase transitions" such as the snowball transition (or its opposite). 
There is evidence that Earth has gone through at least one globally frozen, "snowball" state in 
the last billion years, which it is thought to have exited after several million years because global 
ice-cover shut off the carbonate-silicate cycle, thereby allowing greenhouse gases to build up to 
sufficient concentration to melt the ice. Due to the positive feedback caused by the high albedo of 
snow and ice, susceptibility to falling into snowball states might be a generic feature of water-rich 
planets with the capacity to host life. This paper has two main thrusts. First, we revisit one- 
dimensional energy balance climate models as tools for probing possible climates of exoplanets, 
investigate the dimensional scaling of such models, and introduce a simple algorithm to treat the 
melting of the ice layer on a globally- frozen planet. We show that if a terrestrial planet undergoes 
Milankovitch-like oscillations of eccentricity that are of great enough magnitude, it could melt out 
of a snowball state. Second, we examine the kinds of variations of eccentricity that a terrestrial 
planet might experience due to the gravitational influence of a giant companion. We show that 
a giant planet on a sufficiently eccentric orbit can excite extreme eccentricity oscillations in the 
orbit of a habitable terrestrial planet. More generally, these two results demonstrate that the 
longterm habitability (and astronomical observables) of a terrestrial planet can depend on the 
detailed architecture of the planetary system in which it resides. 

Subject headings: astrobiology - planetary systems - radiative transfer 



X 



1. INTRODUCTION 

Even very mild astronomical forcings can have 
striking influences on the Earth's climate. Although 
the orbital eccentricity varies between ~0 and only 
~0.06, and the axial tilt, or obliquity, between 22.1° 
and 24.5°, these slight periodic changes are sufficient 
to help drive the Earth into ice ages at regular inter- 
vals. Milankovitch (1941) recognized and articulated 
this possibility in his astronomical theory of the ice 
ages. 1 Specifically, Milankovitch posited a causal con- 
nection between three astronomical cycles (precession 
- 23 kyr period, and variation of both obliquity and ec- 
centricity - 41 kyr and 100 kyr periods, respectively) 
and the onset of glaciation/ deglaciation. Though much 
remains to be discovered about the Milankovitch cy- 
cles, they are now generally acknowledged to have been 

Electronic address: dsp@astro.princeton.edu 

1 Prior to Milankovitch, Adhemar (1842) and Croll (1875) 
argued that astronomical forcing influences glaciation. 



the dominant factor governing the ice ages of the last 
several million years (Bcrgcr 1975; Hays et al. 1976; 
Berger 1976, 1978; Berger et al. 2005). 

The nonzero (but, at just 0.04, nearly zero) eccen- 
tricity of Jupiter's orbit is a significant driver of the 
Earth's eccentricity Milankovitch cycle (though the ef- 
fects of the other planets in our solar system are im- 
portant, as well). If Jupiter's eccentricity were much 
greater, it could drive larger amplitude variations of 
the Earth's eccentricity. This same mechanism might 
be operating in other solar systems. Among the more 
than 450 currently known extrasolar planets, there are 
many that have masses comparable to Jupiter's and 
that are on highly eccentric orbits; ^20% of the known 
exoplanets have eccentricities greater than 0.4, includ- 
ing such extreme values as 0.93 and 0.97 (HD 80606: 
e = 0.93, Naef et al. 2001; Gillon 2009; HD 20782b: 
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e = 0.97, O'Toole et al. 2009). 2 Furthermore, tan- 
talizing evidence suggests that lower mass terrestrial 
planets might be even more numerous than the giant 
planets that are easier to detect (O'Toole et al. 2009; 
Mayor et al. 2009). Therefore, it seems highly likely 
that many terrestrial planets in our galaxy experience 
exaggerated versions of the Earth's eccentricity Mi- 
lankovitch cycle. Furthermore, although the Earth's 
obliquity appears always to be contained within a 2.4° 
window 3 , the obliquity of planets that lack large moons 
might vary cyclicly or chaotically, spanning a much 
larger range of angles (Laskar & Robutel 1993; Laskar 
et al. 1993; Neron de Surgy & Laskar 1997). Further- 
more, if Jupiter were significantly closer to the Earth 
(for example, if it were at 3 AU), the Earth's obliq- 
uity would vary chaotically even with the Moon present 
(Williams 1998a, b). In addition, although they are ob- 
servationally unconstrained, models suggest that the 
spin-rates and initial obliquities of exoplanets might 
assume a wide range of values (Agnor et al. 1999; 
Chambers 2001; Kokubo & Ida 2007; Miguel & Brunini 
2010). 

These kinds of changes in external forcing could 
have a dramatic effect on life that requires liquid water. 
Since seminal work by Dole (1964) and Hart (1979), 
a variety of theoretical investigations have examined 
the possible climatic habitability of terrestrial exoplan- 
ets. Kasting et al. (1993) emphasized that the habit- 
ability of an exoplanet depends on the properties of 
the host star. Several authors have considered how a 
planet's climate and climatic habitability depend on 
the properties of the planet, as well (Williams & Kast- 
ing 1997; Franck et al. 2000a,b; Williams & Pollard 
2002; Gaidos et al. 2005; von Bloh et al. 2007; Selsis 
et al. 2007; Dobrovolskis 2007, 2009; von Bloh et al. 
2008; Spiegel et al. 2008, 2009). In particular, two 
recent works have focused on the climatic effect of or- 
bital eccentricity. Williams & Pollard (2003) used a 
general circulation climate model to address the ques- 
tion of how the Earth's climate would be affected by a 
more eccentric orbit. A companion paper to this one 
(Dressing et al. 2010, hereafter D10) uses an energy 
balance climate model to explore the combined influ- 
ences of eccentricity and obliquity on the climates of 
terrestrial exoplanets with generic surface geography. 
A more eccentric orbit both accentuates the ratio of 
stellar irradiation at periastron to that at apoastron, 
and increases the annually averaged irradiation (in pro- 
portion to (1 — e 2 ) -1 / 2 , where e is eccentricity). Thus, 
periodic oscillations of eccentricity will cause concomi- 
tant oscillations of both the degree of seasonal extremes 
and of the total amount of starlight incident on the 
planet in each annual cycle. Since these oscillations 
depend on gravitational perturbations from other com- 
panion objects, the present paper can be thought of 
as examining how a terrestrial planet's climatic habit- 

2 See http: / /exoplanet. eu 

3 Even these small variations can cause significant climatic 
changes (Drysdalc et al. 2009). 



ability depends not just on its star, not just on its own 
intrinsic properties, but also on the properties of the 
planetary system in which it resides. 

There is evidence that, at some point in the last 
billion years, Earth went through a "Snowball Earth" 
state in which it was fully (or almost fully) covered with 
snow and ice (see the review by Hoffman & Schrag 2002 
and references therein; Williams 1993 and Williams 
et al. 1998 suggest an alternative interpretation of the 
data). The high albedo of ice gives rise to a positive 
feedback loop in which decreasing surface temperatures 
lead to greater ice-cover and therefore to further net 
cooling. As a result, the existence of a low-temperature 
equilibrium climate might be a generic feature of water- 
rich terrestrial planets, and such planets might have a 
tendency to enter snowball states (Tajika 2008). The 
ice-albedo feedback makes it quite difficult for a planet 
to recover from such a state (Pierrehumbert 2005). 
In temperate conditions, the Earth's carbonate-silicate 
weathering cycle acts as a "chemical thermostat" that 
tends to prevent surface temperatures from straying 
too far from the freezing point of water (Walker et al. 
1981; Zeebe & Caldeira 2008). A snowball state would 
interrupt this cycle. The standard explanation of how 
the Earth might have exited its snowball state is that 
this interruption of the weathering cycle would have al- 
lowed CO2 to build up to concentrations approaching 
~1 bar over 10 6 -10 7 years, at which point the green- 
house effect would have been sufficient to melt the 
ice-cover and restore temperate conditions (Hoffman 
& Schrag 2002; Caldeira & Kasting 1992). However, 
an exoplanet in a snowball state, that is undergoing 
a large excitation of its eccentricity, might be able to 
melt out of its globally frozen state in significantly less 
time, depending on the magnitude of the eccentricity 
variations and on other properties of the planet. Ex- 
ploring this possibility is a primary focus of this paper. 

But is it reasonable to expect Earth-like planets 
around other stars to have large eccentricities? The 
very circular and well-behaved orbits in the inner So- 
lar System are thought to be the result of dissipative 
processes during the late stages of planetary growth, 
in particular dynamical friction resulting from interac- 
tions between growing protoplanets and remnant plan- 
etesimals (O'Brien et al. 2006; Morishima et al. 2008; 
Raymond et al. 2006, 2009). After the relatively short 
(~100 Myr) phase of planet formation, the Solar Sys- 
tem's long-term evolution is thought to have been dom- 
inated by dynamical interactions between the planets. 
Currently, the Earth's eccentricity oscillates between 
almost zero and about 0.06 with a ^100,000 year peri- 
odicity (Laskar 1988; Quinn et al. 1991). The dynam- 
ics of the Earth's orbit is controlled by secular forcing 
from the other planets and undergoes chaotic evolution 
on long timescales (Laskar 1989, 1990). In the context 
of the known extra-solar planets and a more general 
picture of planet formation and evolution, we expect 
a much wider range of outcomes than what is seen in 
the Solar System (Kita et al. 2010). 
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During the first few million years of planetary 
growth, gaseous protoplanetary disks play an impor- 
tant role in the dynamical evolution of accreting plan- 
ets (see review by Papaloizou & Terquem 2006). Or- 
bital eccentricities of protoplanets may be increased 
by turbulent forcing, interactions between protoplanets 
or with a slightly elliptical gas disk (Papaloizou et al. 
2001; Kokubo & Ida 2002; Oishi et al. 2007; D'Angelo 
et al. 2006) and decreased by a combination of tidal 
damping and the effects of dynamical friction and aero- 
dynamic gas drag (Adachi et al. 1976; Tanaka & Ward 
2004; Kominami & Ida 2004; O'Brien et al. 2006). The 
final eccentricities of a particular system of terrestrial 
planets results from a competition between these nu- 
merous effects. In addition, given that the lifetime of 
gaseous disks is far shorter than the Earth's measured 
accretion timescale (Haisch et al. 2001; Pascucci et al. 
2006; Touboul et al. 2007; Allegre et al. 2008), the or- 
bital configuration of any giant planets in the system 
will play an important role in the final sculpting of the 
terrestrial planets (e.g., Levison & Agnor 2003; Ray- 
mond et al. 2004). 

Gap-opening planets (M p J> Mg aturn ; Crida et al. 
2006) may acquire eccentricities on the order of 0.1- 
0.2 from planet-disk interactions (Goldreich & Sari 
2003; D'Angelo et al. 2006). However, these values 
are too small to explain the eccentricity distribution of 
the known extra-solar giant planets, which is centered 
at 0.2-0.3 and contains values higher than 0.9 (Butler 
et al. 2006; Wright et al. 2009). Several mechanisms 
have been proposed to explain the observed eccentric- 
ity distribution (see Ford & Rasio 2008 for a summary). 
Currently, the leading candidate is the planet-planet 
scattering mechanism, which proposes that most exo- 
planet systems (75-100% of them) became dynamically 
unstable after their formation, leading to close encoun- 
ters between planets and the eventual destruction of 
some planets (Rasio & Ford 1996; Weidenschilling & 
Marzari 1996; Juric & Tremaine 2008; Thommes et al. 
2008). Planet-planet interactions between Earth-sized 
objects are much more likely to result in eccentricity- 
damping collisions than eccentricity-increasing scatter- 
ing events (Goldreich et al. 2004). However, scatter- 
ing among giant planets in the vicinity of rocky pro- 
toplanets can destroy the building blocks of Earth-like 
planets if the giant planets are close-by, and the pro- 
toplanets that survive tend to have large eccentrici- 
ties (Veras & Armitage 2006). In addition, terrestrial 
planets that form in systems with very massive or very 
eccentric giant planets tend to themselves have signifi- 
cant eccentricities (Chambers & Cassen 2002; Levison 
& Agnor 2003; Raymond et al. 2004). Simulations of 
terrestrial planet formation show significant variations 
between systems with similar initial conditions (e.g., 
Raymond et al. 2004; Quintana et al. 2007). Given 
that the range of reasonable input parameters (giant 
planet orbits, disk properties) is far wider than that 
tested with simulations, it seems reasonable to expect 
a large diversity in the orbits of extrasolar Earth-like 



planets. 

The rest of the paper is organized as follows: In 
§2, we present our energy balance climate model and a 
dimensional analysis of its behavior. In §3, we discuss 
some scenarios in which our model suggests that there 
might be climatic consequences to eccentricity oscilla- 
tions that would be interesting from the perspective of 
habitability. In §4, we use secular perturbation theory 
and an N-body integrator to investigate what sort of 
planetary system architecture might lead to large am- 
plitude variations in a terrestrial planet's eccentricity. 
Finally, in §5, we conclude. 

2. CLIMATE MODEL 

We use nearly the same 1-dimensional time- 
dependent energy balance model (EBM) that wc have 
used in previous studies of exoplanet climatic habit- 
ability (Spiegel et al. 2008, hereafter SMS08; Spiegel 
et al. 2009, hereafter SMS09; see references therein for 
a justification of this model in habitability studies). 
The present work's model has a slight modification, de- 
scribed in §2.2, to handle starting in a snowball state. 
Our model treats the redistribution of heat as a diffu- 
sive process, and it (or close variants) has been used 
by other authors both in the context of exoplanet hab- 
itability (Williams & Kasting 1997, hereafter WK97) 
and in other contexts, including studies of both Mar- 
tian climate (Nakamura & Tajika 2002, 2003) and the 
Earth's climate (North et al. 1981). 4 The response of 
this model to nonzero eccentricity is explored exten- 
sively in a companion paper (D10). Here, we provide 
a brief synopsis of the model. 

In accord with WK97 and our own previous work, 
we use the following equation for planet surface tem- 
peratures (T) as a function of time (t) and location 
(*): 

C f -£0* "*">£) -*<>-*>-'. « 

where x = sin A and A is latitude. S[x,t] is the diur- 
nally averaged stellar irradiation at the latitude band 
represented by x, on the date represented by t. The 
values of heat capacity C and effective diffusivity D, 
and of the albedo and infrared cooling functions A 
and / are as described in SMS08 and SMS09. We 
briefly recapitulate. The heat capacity has the follow- 
ing values over land (C;), ice (Ci), and ocean (C D ): 
Q = 5.25 x 10 9 erg cm" 2 K"\ C t = 9.2Q when 263 < 
T < 273 K and 2.0C, when T < 263 K, and C a = 40C/; 
these values are summarized in Table 1. The diffusivity 
D = 5.394 x 10 2 erg cm~ 2 s" 1 K" 1 . The albedo is A 2 
of SMS08 and SMS09: A[T] = 0.525 - 0.245 tanh[(T - 
268 K)/5 K]. We adopt two different functional forms 
for the infrared cooling. The first is I 2 of SMS08 and 
SMS09: I 2 [T] = aT 4 / {l + 0.5925(T/273 K) 3 }, where 

4 Suarcz & Held (1979) used a more sophisticated analog of 
this EBM, but in an Earth-centric context. 
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a = 5.67 x 1(T 5 crm cirT 2 s" 1 K~ 4 is the Stefan- 
Boltzmann constant. The second is Iwk97 of SMS09 
(and WK97), which depends not only on T but on 
Pco 2 : the CO2 partial pressure. This function, pre- 
sented in the appendix of WK97, is a parameterized 
fit (third order in T and fourth order in ln[pco 2 ]) to 
results from a radiative-convective model (Kasting & 
Ackerman 1986; Kasting 1988, 1991). We solve equa- 
tion (1) on a grid of 145 points equally spaced in lat- 
itude. We use a time-implicit numerical scheme and 
an adaptive time-step, as described in SMS08 and in 
Hameury et al. (1998). Obliquity set to Earth-like 
23.5° and to 90°, and initial temperature is typically 
set to be 100 K, with the melting of ice handled as 
described in §2.2. 

2.1. Dimensional Analysis of EBM 

As written in equation (1), the energy balance equa- 
tion appears similar to a diffusion equation, but is writ- 
ten in terms of a nondimensional spatial variable (x). 
It is instructive to rewrite equation (1) both in fully 
dimensional and in fully nondimensional form. 

The (1-D, latitudinal) forced heat equation on the 
sphere may be written as 

7/ — T TT" S cos \s I R] — — } = — \S(l - A) - 1} . 

dt cos[s/R] ds\ 1 ' J ds J C 1 v ; s 

Here, R is the planetary radius, s = RX is the distance 
north of the equator, and K yy is a dimensionally proper 
meridional eddy diffusion coefficient. This diffussivity 
may be expressed as K yy = R 2 /rdig, where Tdis is the 
diffusion time. The factors of cos[s/R] in equation (2) 
are metric terms that arise from writing the diffusion 
equation in spherical coordinates. 

Equation (2) may be converted to equation (1) by 
replacing s with x = sin[s/i?], multiplying by C, and 
defining D = C /r^is- For the values of C and D quoted 
above (WK97, SMS08), the diffusion timescale may be 
written as 

r di ffw(C/C,)(9.7x 10 6 s) 

w(C/Ci) (4 months). (3) 

We find that over land, thick ice, thin ice, and ocean, 
Tdiff has the following respective values: ^4 months, 
^7 months, ~3 years, and ~12 years. The thermal 
diffusivities (K yy ) in Table 1 (<~ 10 9 cm 2 s _1 for atmo- 
sphere above ocean, 40 times greater for atmosphere 
above land) might seem surprisingly large for Earth; 
in fact, however, they comport with what is expected 
when a process that involves large-scale advective mo- 
tions is modeled as being purely diffusive (Kao & Kau 
1980; Keeling & Heimann 1986). 5 

5 Lorenz (1979) justified the diffusive approximation for stud- 
ies of large scale climatic trends; this treatment has been used 
in many EBMs in the geophysical literature in the last several 
decades (e.g., Suarez & Held 1979; North et al. 1981, 1983). 
Furthermore, Showman et al. (2009) provide an illuminating dis- 



Alternatively, we may write a fully nondimensional 
version of equation (1) by mapping t 1— > P rb^*, T n- 
T T*, J H J r, and S ^ I S*, with I = <tT 4 (ct 
is the Stefan-Boltzmann Constant). Here, P rb is the 
orbital period and t* is a dimensionless time variable; 
T is a typical temperature and T* is a dimension- 
less temperature variable; Iq is the blackbody cool- 
ing rate associated with temperature To and I* is a 
dimesionlcss infrared cooling variable; and S* is a di- 
mensionless insolation variable. Noting that the ra- 
diative timescale, or thermal inertia, may be written 
r rad = CT /I = (C/a)T ~ 3 , we are left with the fol- 
lowing nondimensional form of the equation: 




= n*{s*(i-A)-i*} . 



(4) 

Here, JC* = P rb/7"diff may be thought of as a dimen- 
sionless thermal diffusivity, while 1Z* = Porb/^rad may 
be thought of as a nondimensional "thermal elastic- 
ity." If 1Z* 3> 1 the planet's temperature will tend to 
be determined by the instantaneous irradiation; con- 
versely, if 72.* <C 1, the temperature will change little 
over the course of an annual cycle. Writing the equa- 
tion as in (4) makes clear that increasing the orbital 
period (i.e., increasing the semimajor axis a or decreas- 
ing the stellar mass at fixed luminosity) increases the 
model planet's ability both to redistribute thermal en- 
ergy in a given fraction of a year, and to approach ra- 
diative equilibrium. While this is intuitively obvious, it 
is important to notice that, even while holding many of 
a model planet's dimensional parameters constant (in 
particular, while holding C and D fixed), if the orbital 
separation changes then the nondimensional diffusiv- 
ity and elasticity change as well. This is another way 
of looking at an issue that was explored in SMS08: in 
short, timescales matter, and a closer planet is not sim- 
ply a more strongly irradiated version of a more distant 
planet. 

2.2. Modeling a Cold Start 

We place a frozen planet in a variety of pre-set or- 
bits in order to explore the capacity for the orbital 
configuration (specifically, the semimajor axis, the ec- 
centricity, and the obliquity) to thaw such a planet. 
In previous work using this model, we have set "hot 
start" initial conditions in which the initial tempera- 
ture is set far above the freezing point of water. In 
modeling the melting of a Snowball Earth planet, we 
are in a different regime. We have previously implic- 
itly assumed that the latent heat involved in melting 
ice and freezing water is negligible. Is this assumption 
still valid in the case of a snowball planet that might 
have a layer of ice that is a kilometer or more thick? 

The latent heat of melting ice (L- lcc = 3.3 x 
10 9 erg g -1 ) may be considered negligible if it is small 

cussion of the applicability of the diffusive approximation. We 
explored the consequences of varying D in SMS08, SMS09, and 
D10. 
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compared with the mean specific energy that is de- 
posited into ice in a year, (dE/dm). For a crude upper 
bound on this quantity, we consider just incident radi- 
ant energy: 

/ dE\ forb-Frad 

\dm/ hp ice 

where F la( i is the typical incoming radiative flux on a 
planet around a Sun-like star, ft is the height of a layer 
of ice, and p- lcc is the density of ice. At a distance a 
around a Sun- like star, P orb = (1 yr)(a/l AU) 3/2 , and 
the typical stellar irradiation flux is 



^~* /4 (nu)" 



:3.4 x 10 5 



(l Au) 



erg cm 2 s 1 , 



(6) 



where qo is the solar constant, the flux at the subsolar 
point on Earth. Therefore, L lcc (dE/dm) when 



h Z 



^orb -^rad 
Picc-^icc 



-1/2 
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Since the thickness of ice on Earth varies seasonally by 
far less than 35 m (Semtner 1976), it is not surpris- 
ing that SMS08 and SMS09 found excellent agreement 
between our EBM and the Earth's seasonal variations 
of temperature (SMS08) and radiative and diffusive 
fluxes (SMS09). A planet at 1 AU from a Sun- like 
star, however, does not receive enough radiant energy 
in a year to melt a layer of ice that is more than 35 m 
thick. 

We here adopt a simple algorithm to handle the 
melting of an ice layer. 6 At a given location, if the 
temperature has never reached 273 K, the albedo is 
set to its maximum value (0.77), to indicate that the 
surface is still covered with snow and ice. When the 
local temperature reaches 273 K, further excess of ra- 
diant / diffusive energy (per area) of AEa is treated as 
changing the thickness of the local ice layer by Aft,, 
where Aft = —AEA/(pi C eLi C e)- During such periods 
of melting, the albedo reduces to 80% of its maximum 
value (A = 0.616) because snow and ice are less shiny 
when melting than when fully frozen (K0ltzow 2007; 
Mellor & Kantha 1989; Curry et al. 2001). Consider 
a patch of surface that has just been melting but, be- 
cause of the onset of winter, now receives less incident 
flux than its radiant and diffusive losses. First, the 
summer's melt-water re-freezes (AEa is negative, so 
Aft, is positive) while the temperature remains fixed at 
273 K. Once all melt-water has re-frozen, subsequent 
energy deficits result in the temperature decreasing, as 
specified by equation (1). 

Instead of assuming an initial ice thickness and in- 
tegrating until all of the ice has melted, we recognize 

6 See Semtner (1976) for a more sophisticated model of the 
vertical diffusion of heat through an ice layer. 



that if a portion of a model planet, during a given year, 
receives more incident radiant energy than its radiant 
and diffusive losses, then it will continue to receive an 
excess of energy in all following years. Therefore, we 
define a patch of surface as having "melted" once there 
has been positive "net melting" over the course of a 
year (i.e., the summer's melt- water does not fully re- 
freeze during the winter). In this way, we compress 
what would be perhaps ^10 3 years of melting on a real 
planet to a single year. If an initially frozen model 
planet's surface has partially melted within the inte- 
gration time of 130 model years, we note this model's 
orbital configuration as one that could melt a planet 
out of a snowball state. Our simple algorithm can be 
criticized for a variety of reasons (e.g., perhaps the 
thickness of the ice ought to be constrained to be con- 
tinuous or differentiable) but, given other uncertainties 
involved in modeling the climate of exoplanets, ours is 
a reasonable first ansatz to explore the principle under 
consideration in this paper. 8 

Finally, over the ^10 5 years that an exo-Earth's 
eccentricity Milankovitch cycle might last, a globally 
frozen planet that is comparably geologically active to 
the Earth might release ~0.01 bars of CO2 due to vol- 
canism (WK97; Holland 1978). With this in mind, we 
examine our model's behavior with three different in- 
frared cooling functions: (i) I 2 of SMS08 and SMS09, 
which is for an Earth-like greenhouse gas abundance 
(-3 x 10~ 4 bars C0 2 ); (ii) 7 W K97 of WK97, with the 
partial pressure of CO2 set to 0.01 bars; and (iii) iwK97 
with variable CO2 partial pressure. Once a cold-start 
model's ice-cover has completely melted somewhere, 
silicate weathering might recommence, thereby remov- 
ing CO2 from the atmosphere. On a real planet, the 
rate at which melting continues would involve a compli- 
cated interplay between the thickness of the ice cover, 
the rates of both CO2 sequestration (through weath- 
ering) and release (through volcanism), and the rate 
of change of astronomical forcing. Since none of these 
values are constrained for generic exoplanets, we in- 
stead adopt a crude version of a "chemical thermostat" 
for cooling function (iii): p C o 2 = 10~ 2 ~( T ~ 250 )/ 27 bars 
if 250 K < T < 290 K; p C o 2 = 3.3 x 10~ 4 bars if 
T > 290 K; and p C o 2 = 0.01 bars if T < 250 K. 9 
In models in which the CO2 concentration is adjusted 
with temperature, this adjustment begins after the the 
temperature has somewhere exceeded 273 K for an en- 
tire model year. 

7 The estimate of ~10 3 years comes from estimating the 
length of time required to melt a <~kilometer-thick ice layer at 
~100 cm yr _1 , though cither the thickness or the rate of melting 
could differ significantly from these estimates. 

8 See WK97 for an extensive discussion of the utility and 
limitations of energy balance climate models in the context of 
exoplanet climate modeling. 

9 The work of Le Hir et al. (2009) suggests that it might 
take significantly longer than we assume for weathering processes 
to reestablish equilibrium CO2 concentrations after a snowball 
event. 
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3. CLIMATE RESULTS 

When a terrestrial planet's orbital eccentricity or 
obliquity varies with time, a variety of climate effects 
are possible. The climate simulations described in Df 0, 
Williams & Pollard (2002), and Williams & Pollard 
(2003) give us a rough idea of how the climates of some 
planets with time- variable astronomical forcing might 
behave. But if a planet undergoes a phase transition 
- entering into a snowball state, or melting out of one, 
the effect of varying eccentricity is less clear. 

Here, we address two questions, (i) What eccentric- 
ity is required to melt out of a snowball state? (ii) How 
does an oscillating eccentricity (a Milankovitch-like cy- 
cle) affect the climate of a planet that has slipped into 
a snowball state? 

3.1. Cold-Start Planets 

Rather than attempting to map fully the region of 
parameter space that leads to melting planets, we in- 
stead explore the eccentricity required to melt a cold- 
start model in several illustrative examples. In partic- 
ular, we consider snowball planets at 0.8 AU and at 
1.0 AU. 

For planets with nonzero eccentricity, two angles 
are required to describe the direction of the spin-axis 
with respect to the orbital ellipse. D10 defines these 
as two obliquity angles: the "polar obliquity" (0 P ), 
which is the familiar angle between the spin axis and 
the orbital plane's normal, and the "azimuthal obliq- 
uity" (0 a ), which is the difference in true anomaly from 
the periastron point to the northern winter solstice. 10 
D10 suggests that for planets with high eccentricity, 
azimuthal obliquity angles near 30° might be some- 
what more stable against global freeze-over. Taking 
our cue from this result, we set the azimuthal obliq- 
uity angle to 30° for all models considered here. In our 
brief survey of parameter space, we consider two val- 
ues of polar obliquity: Earth-like 23.5°, and an extreme 
value of 90°. We note that, if an exoplanet experiences 
Milankovitch-like variations of precession and obliquity 
(and in particular if its polar obliquity is not stabilized 
by a moon), then its azimuthal and polar obliquity an- 
gles might be expected to swing through a large portion 
of all possible values. 

The eccentricities required to melt out of a global 
snowball state are shown for each of 8 model planets 
in Table 2 (two semimajor axes x two polar obliquities 
x two cooling functions). Not surprisingly, the eccen- 
tricity need not reach values that are as great if the 
greenhouse effect is enhanced by greater atmopsheric 
CO2 content. The values listed (ranging between 0.29 
and 0.86) may all be excited by Milankovitch-like cy- 
cles in other solar systems (see Figs. 3 and 4). If the 
geochemistry of an exoplanet should cause either its 

1(1 Equivalcntly, it is the angle between the vector from the 
star to the periastron point and the projection of the planet's 
spin-vector onto the orbital plane. 



ice to have lower albedo than we assumed, or its vol- 
canos to release CO2 significantly more rapidly than 
we considered, the required eccentricity values might 
be even lower. Furthermore, it is likely that a real 
planet would have additional feedback processes that 
might either amplify or damp the warming caused by 
increased eccentricity. 

Figure 1 provides a graphical demonstration of the 
temperature evolution on cold-start planets. Both the 
left and right panels depict the climate evolution of 
models with polar obliquity of 23.5°, azimuthal obliq- 
uity of 30°, and eccentricity of 0.8. Both models use 
the WK97 infrared cooling function. In the left panel, 
the CO2 partial pressure is held constant at 0.01 bars. 
In the right panel, the CO2 partial pressure begins at 
0.01 bars but then, a year after there is positive net 
melting somewhere on the planet, it adjusts with tem- 
perature as described in §2.2. For the first ~9 years, 
both panels evolve identically. After ~9 years, an equa- 
torial region of year-round temperate conditions devel- 
ops. The model in the left panel at this point contin- 
ues to warm progressively. Once the equatorial region 
melts, the fraction of surface that has melted ice-cover 
grows steadily until the entire planet has melted, and 
temperatures eventually grow to more than 400 K over 
much of the planet (not shown). In the right panel, in 
contrast, the thermostat begins to operate after ^9 
years, causing the CO2 concentration to decrease sig- 
nificantly. The negative feedback causes this model to 
achieve a stable equilibrium climate. 

We demonstrate in §4 that, depending on system 
architecture, large variations of eccentricity are possi- 
ble on timescales of less than 10 6 years. In this con- 
text, Table 2 and Fig. 1 make it clear that there are 
possible circumstances in which eccentricity excitation 
could dramatically warm an initially frozen planet. In- 
triguingly, in order for a planet to melt out of a snow- 
ball state by virtue of increased eccentricity and then 
to remain with a stable, habitable climate, it might 
have to have released some greenhouse gas during its 
frozen period. In other words, the classical explanation 
for how the Earth might have exited its own snowball 
state might need to apply to some extent, although 
the increased eccentricity can significantly reduce how 
much greenhouse gas is required. 

3.2. Experimenting with a Pseudo-Milankovitch Cycle 

On a real planet with an extreme Milankovitch- 
like cycle of eccentricity oscillations, the climate never 
quite reaches a stable year-to-year equilibrium, because 
the eccentricity is constantly changing. Still, even for 
the most extreme cases explored in §4, the changes 
in eccentricity are very small from one year to the 
next, since the periods of eccentricity oscillation are 
?Z 10 3 years. As a numerical experiment to probe what 
can happen when eccentricity is varied as a climate 
model runs, we compress two Milankovitch-like eccen- 
tricity cycles into 25 years. One cycle is for a planet 
at 1 AU, and the other is for a planet at 0.8 AU. For 
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simplicity, we set the eccentricities of both model plan- 
ets to vary sinusoidally. The 1-AU planet has an ex- 
treme Milankovitch-like cycle, with eccentricity vary- 
ing from to 0.83: e = 0.83{1 + 0.5 sin[(i- 10)2tt/25]}, 
and the 0.8-AU planet has a more typical cycle, with 
eccentricity varying from 0.1 to 0.33. We start both 
planets with temperate conditions, warm equator and 
cold poles (T init [A] = 240 + 60 cos 2 A). Furthermore, 
when the eccentricity is below 0.05 for the 1-AU planet 
and below 0.125 for the 0.8-AU planet, we set the 
albedo to 0.8, simulating a catastrophic event that 
plunges the planet into a snowball state in which in- 
frared cooling is governed by JwK97, initially with con- 
stant pco 2 — 0-01 bars. Figure 2 displays the results of 
four such model runs (two for planets at 1 AU, two for 
planets at 0.8 AU; all planets have 23.5° polar obliq- 
uity). In the left panels, the CO2 concentration is held 
constant at 0.01 bars throughout; in the right panels, 
the CO2 concentration drops according to the prescrip- 
tion in §2.2 once a region of surface has remained above 
freezing for an entire year. As in Fig. 1, the negative 
feedback of the thermostat's turning on mutes the pos- 
itive feedback of the ice's melting, resulting in a smaller 
fraction of the ice-cover melting within the cycle. 

3.3. Discussion 

In §3.1 and §3.2, we presented examples of plan- 
ets with eccentricities between 0.1 and 0.86, and with 
eccentricity oscillation cycles of magnitude 0.83 and 
0.23. Although oscillations greater than 0.5 in a ter- 
restrial planet's eccentricity (such as those in Fig. 2 and 
some of those in Table 2) might be considered extreme, 
oscillations of ^0.3 may be routine. Furthermore, a 
planet's variable eccentricity need not always return to 
zero (see the lower right panel of Fig. 3), so our galaxy 
might host many planets whose eccentricities vary be- 
tween, say, 0.1 and 0.4. Given this, it is worth keeping 
in mind that the snowball cases considered above might 
be extreme, as well. Some climate simulations suggest 
that the Earth's snowball state(s) (if in fact it went 
through any) might have left an equatorial band of 
open water, perhaps because intense winds (driven by 
the sharp temperature contrast between frozen and un- 
frozen tropical ocean) stirred up the typically ^50 m 
wind-mixed layer (Hartmann 1994) to much deeper, 
tapping into the large heat capacity of the deep ocean 
(Hyde et al. 2000). The idea of a temperate region 
persisting throughout the (partial) snowball state is 
also supported by the fact that life evidently survived 
any such states, which might have been difficult if the 
whole surface was covered with ice for ^ 10 6 years. 
Much milder (and more probable) eccentricity excur- 
sions could cause a planet in either a partial snowball 
state or in an analog of the Earth's ice ages to warm 
significantly. 

In our models, we have held each planet's obliquity 
constant (at 23.5° or 90°, depending on the model). 
This assumption is purely to limit the complexity of 
our modeling endeavor; it is not motivated by astro- 



physical intuition or theory. The architecture of a plan- 
etary system affects not only a terrestrial planet's or- 
bital eccentricity, as considered in §4, but also its obliq- 
uity. Since terrestrial exoplanets' obliquities might 
vary significantly (Williams 1998a, b), our assumption 
of constant obliquity restricts the richness of the prob- 
lem; the evolution of climate and ice-distribution on a 
planet whose spin-axis direction and orbital eccentric- 
ity are both varying wildly is much more complex than 
we have considered. 

Terrestrial exoplanets might someday be discovered 
that have spectroscopic and photometric character- 
istics indicative of being covered with snow and ice 
(and possibly of being shrouded in cold CO2 clouds) - 
i.e., that seem to be in snowball states. If such plan- 
ets are found, analyzing the full architectures of their 
planetary systems will allow us to understand whether 
their eccentricities might vary enough (and on what 
timescales) to help them to melt out of their frozen 
states. It might even be possible to determine whether 
a planet is in the process of increasing (or decreasing) 
its orbital eccentricity, and therefore whether it could 
be in the process of melting (or freezing), which might, 
at some point, be testable through direct observations. 

4. ORIGIN OF EARTH-LIKE EXOPLANET 
ECCENTRICITIES 

For the purposes of this study, we are interested 
in the eccentricity evolution of Earth-like planets, in 
terms of both the amplitude of oscillation and the pe- 
riod of oscillation. Given the importance of resonances 
and the potential for chaos, it is impossible to under- 
stand the detailed dynamical evolution of a planetary 
system without full knowledge of the system's archi- 
tecture. However, we can make simple assumptions to 
get an idea of the range of possibilities. We therefore 
restrict ourselves to the simple case of an Earth analog 
at 1 AU evolving under the perturbations from a sin- 
gle massive body, in most cases a Jupiter-mass giant 
planet. To study these simple, two-planet systems we 
numerically integrated the orbits of several systems to 
look at the range of outcomes and correlations between 
the perturbing body's mass/orbit and the evolution of 
the Earth-like planet. 

We used two different numerical techniques to in- 
tegrate our Earth analog - single perturber systems. 
First, to explore the range of outcomes in coplanar sys- 
tems, we integrate the orbit-averaged rates of change 
for the planets' eccentricities e and longitudes of peri- 
center zn using the equations of Mardling & Lin (2002) 
and Mardling (2007). Our second numerical tech- 
nique was to directly track the 3-dimensional orbits 
of the planets with the Mercury symplectic integrator 
(Chambers 1999). We used the orbit-averaged method 
to study coplanar systems, and Mercury to verify the 
results of the orbit-averaged code and also to study 
mutually inclined systems. We evolved a variety of 
systems for 1 Myr, in each case including an Earth- 
mass planet at 1 AU (initially on a circular orbit), and 
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neglecting the effects of general-relativistic precession. 
Our simulated systems were chosen not to make predic- 
tions about specific known cxoplanetary systems but 
rather to illustrate the variety of eccentricity oscilla- 
tions of terrestrial planets in dynamically stable (al- 
beit simplified) planetary systems, and to show simple 
correlations between the properties of the perturbing 
body and the eccentricity of the Earth analog. 

The results of our integrations are summarized in 
Fig. 3. The first three panels show the effect of a partic- 
ular system parameter: the giant planet mass (top left 
panel), the giant planet eccentricity (top right panel) 
and the giant planet's orbital distance (bottom right 
panel). The amplitude of the Earth-like planet's ec- 
centricity oscillation is a function of the giant planet's 
semimajor axis and eccentricity but is independent of 
the giant planet's mass (top left panel). However, the 
oscillation frequency (often referred to as the secular 
frequency) increases linearly with giant planet mass. 
The amplitude of the Earth analog's eccentricity os- 
cillation decreases for lower giant planet eccentricities, 
and the period of oscillation increases with lower giant 
planet eccentricities (top right panel of Fig. 3). 

The bottom left panel of Fig. 3 shows two systems 
that both induce eccentricity oscillations of about 0.1 
in the Earth analog, but on very different timescales 
of 3,000 and 170,000 years. The two systems are 
markedly different in their architecture: the fast- 
oscillating case contains a giant planet interior to the 
Earth- like planet (giant planet at 0.5 AU with an ec- 
centricity of 0.1); and the slow-oscillating case contains 
an exterior giant planet (at 5 AU with an eccentricity 
of 0.4). 

The first three examples from Fig. 3 were chosen 
such that the giant planets' orbits are consistent with 
the formation of an Earth-like planet at 1 AU in sys- 
tems with both inner (Raymond et al. 2006; Mandcll 
et al. 2007) and outer (Raymond 2006) giant planets. 
For the fourth example, we searched for extreme cases 
that would be dynamically stable, yet induce very large 
amplitude oscillations in the Earth analog's eccentric- 
ity. One potential source of large eccentricities comes 
from a consideration of inclinations called the Kozai 
mechanism (Kozai 1962). For a star-planct-companion 
system, where the companion can be stellar or plane- 
tary, an inclination between the planet and compan- 
ion's orbits / larger than 39.2° can induce eccentricity 
(and inclination) oscillations in the planet's orbit, with 
a maximum eccentricity e max of 

e max ~ \jl - ^cos 2 (7) , (8) 
and with a periodicity Pk oz of roughly 




In equation 9, M* is the stellar mass, a is semimajor 
axes, and subscripts 1 and 2 refer to the planet and the 



companion, respectively (Innanen et al. 1997; Holman 
et al. 1997; Ford et al. 2000). For reasonable estimates 
of the mutual inclination between planetary and com- 
panion orbits, the majority of binary stars should in- 
duce large eccentricity oscillations in any planets that 
form around the primary star (Takeda & Rasio 2005). 

Two extreme cases are shown in the bottom right 
panel of Fig. 3, representing large-scale eccentricity 
forcing from two dynamical sources: secular forcing in 
the coplanar system, as seen in the previous examples, 
and the Kozai mechanism. The secular forcing exam- 
ple is a system with a 10 Jupiter-mass outer planet at 
30 AU with an orbital eccentricity of 0.925. 11 In this 
case, the Earth analog started with an eccentricity of 
0.45; thus, the "backstory" for this system involves 
an impulsive perturbation to the Earth analog that in- 
creased its eccentricity considerably, followed by stable 
secular interactions with the giant planet. The Kozai 
mechanism example contains a 10 Jupiter-mass planet 
orbiting at 10 AU with an eccentric and inclined orbit: 
eccentricity of 0.25 and an inclination of 75° with re- 
spect to the Earth- like planet's starting orbital plane. 
In both cases, the eccentricity of the planet at 1 AU 
reaches 0.9 or higher on a timcscale of one to several 
hundred thousand years. However, it is important to 
note that the inclination variations of these planets are 
markedly different. The system for which the Kozai ef- 
fect is important induces very large variations in the 
Earth analog's inclination, going so far as to sometimes 
rotate in a retrograde sense with respect to its initial 
orbital plane. In contrast, the secularly forced system 
remains virtually coplanar. 

Resonances and the proximity of the giant planet 
are both of great importance to the eccentricity evolu- 
tion of an Earth-like planet. To test these effects we 
performed a small number of integrations of systems 
containing one giant planet and a disk of massless test 
particles whose behavior represents an ensemble of hy- 
pothetical systems with Earth-like planets at different 
orbital distances. We used the hybrid version of the 
Mercury integrator, taking care to choose a step size 
small enough to resolve our innermost orbits at least 
20 times. Figure 4 shows the maximum eccentricity 
reached by test particles in the terrestrial planet zone 
for two different orbital configurations of a Jupiter- 
mass planet, one case with the giant planet exterior 
to 1 AU (shown in black) and one case with the giant 
planet interior to 1 AU (shown in grey). As expected, 
the eccentricity decreases farther from the giant planet, 
in opposite radial directions for the two cases, which 
induce similar eccentricities in planets at 1 AU. The 
"blips" in each curve correspond to individual mean 
motion resonances with the giant planet, which are 
stronger than in the Solar System because of the giant 
planets' eccentricities (e.g., Murray & Dermott 2000). 
The 4:1 resonance with the inner giant planet is visible 

11 Note that three known planets (HD 4113b, HD 80606b, 
and HD 20782b) have orbits with eccentricities of 0.9 or larger 
(Tamuz et al. 2008; Naef et al. 2001; Jones ct al. 2006). 
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at 1 AU. 

To summarize, the eccentricities of Earth-like plan- 
ets are determined cither during their formation or the 
subsequent dynamical evolution of the planetary sys- 
tem - in many cases this will include a dynamical in- 
stability among the giant planets. Once the system 
has settled down, the eccentricity of terrestrial planets 
is determined by the influence of massive perturbers 
such as giant planets and binary companions, via sec- 
ular, resonant, and Kozai forcing. For simple systems 
with a single perturbing planet, the eccentricity os- 
cillation amplitude and period arc determined by the 
giant planet's mass, semimajor axis, and eccentricity. 
For large mutual inclinations, the Kozai mechanism 
is governed by the inclination between the planetary 
and companion orbits, and the mass of the compan- 
ion. The main difference between systems with giant 
planets interior to the terrestrial planets and systems 
with giant planets exterior to the terrestrial planets is 
that the eccentricity oscillation timescale (the secular 
timescale) is much shorter for close-in giant planet sys- 
tems. In systems that do not form any giant planets 
and have no distant companions, there is no obvious 
mechanism by which large eccentricities should be in- 
duced in Earth-like planets. We therefore expect plan- 
ets in such systems to have relatively circular orbits 
(Raymond et al. 2007). 

The so-called dynamical habitability of planetary 
systems has been studied via the stability of test parti- 
cles in the habitable zone of the known extra-solar sys- 
tems (e.g., Menou & Tabachnik 2003; Barnes & Ray- 
mond 2004; Rivera & Haghighipour 2007). This is an 
interesting exercise because it can help observers by 
constraining the orbital location of additional planets 
(see Barnes et al. 2008) and can also test the eccen- 
tricity variations of stable test planets. Such stable 
planets could indeed exist in these systems below the 
detection threshold and represent Earth analogs. 

5. CONCLUSION AND SUMMARY 

We have demonstrated that suitably structured ex- 
oplanetary systems could contain Earth-like planets 
that undergo exaggerated versions of the Earth's eccen- 
tricity Milankovitch cycles. Plausible arrangements of 
giant and terrestrial planets could lead to an Earth-like 
planet's eccentricity oscillating between and more 
than 0.9, on time-scales of a few times 10 3 years to 
a few times 10 5 years. Oscillations such as these could 
have important consequences for climate and habitabil- 
ity. 

A companion paper (D10) investigates the influence 



of a temperate terrestrial planet's eccentricity on its 
habitability. Here, we have examined how eccentric- 
ity could affect a globally frozen terrestrial planet by 
introducing an algorithm to allow our EBM to treat 
the melting of a global ice layer. We found that in- 
creasing the eccentricity of an initially-frozen model to 
sufficiently high values can cause the model to melt out 
of the snowball state. Such high values of eccentricity, 
however, might place a planet at risk of heating up 
to temperatures greater than 400 K (perhaps even un- 
dergoing a runaway greenhouse effect) if there is not a 
concomitant negative feedback, such as from a chemi- 
cal thermostat similar to the Earth's carbonate-siicate 
weathering cycle. Because the climatic influence of in- 
creased eccentricity depends on features of a planet's 
geochemistry that might differ from planet to planet, 
our results should not be taken as the final word on 
precisely what values of eccentricity would be required. 
Nevertheless, the principle remains that if a planet's 
eccentricity is excited to large values, it is conceivable 
there are situations in which this might actually in- 
crease, not decrease, its habitability. 

The longterm habitability of a planet depends on 
the properties of its star, on the properties of its own 
orbit, spin, and geochemistry, and also on the arrange- 
ment of other companion planets. When, in the com- 
ing years, we find Earth-sized planets at orbital sep- 
arations that could be consistent with habitable cli- 
mates, it will be important to study the architectures 
of those systems, so as to determine what kinds of 
Milankovitch-like cycles might govern the longterm cli- 
matic habitability of such worlds. 
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TABLE 1 

Parameter Values Over Land, Ocean, Ice 



Parameter 


Thick Ice 
(J < 263 K) 


Thin Ice 
(263 < i < 273 K) 


Land 


Ocean 




2Ci 


9.2C ; 


Ci = 5.25 x 10 9 erg cm" 2 K" 1 


40C, 


T"diff b 


7 months 


2.8 years 


3.7 months 


12 years 


T ra d C 


3 months 


1.1 years 


1.4 months 


5 years 


K d 
L ^yy 


2.1 x 10 10 cm 2 s" 1 


4.5 x 10 9 cm 2 s" 1 


4.2 x 10 10 cm 2 s" 1 


1.0 x 10 9 cm 2 s" 1 




1.6 


0.3 


3 


0.08 


TV i 


4 


0.9 


8 


0.2 



a C is the heat capacity. 

b Tdiff = C/D is the diffusion time, where D = 5.394 X 10 2 erg cm~ 2 s _1 K _1 . 

c T r ad = (C/o-)T -' A = (1.4 months) (C/Ci)(T/To)" 3 is the radiative time, taking T to be an Earth-like 290 K. 

d K yy = -R© 2 /rdiff is the eddy diffusion coefficient. 

C K* = Porb/ T diff is the nondimcnsional diffusivity of equation (4). 

f 7£* = Porb/^rad is the nondimcnsional thermal elasticity. 



TABLE 2 

Eccentricity Required to Melt Cold-Start Models 



Semimajor Axis 
(AU) 


8 P 

(degrees) 


9a 

(degrees) 


IR Cooling Function 


Required Eccentricity 


0.8 


23.5 


30 




0.39 


0.8 


90 


30 


h 


0.51 


1.0 


23.5 


30 


h 


0.86 


1.0 


90 


30 


h 


0.79 


0.8 
0.8 
1.0 
1.0 


23.5 
90 

23.5 
90 


30 
30 
30 
30 


/WK97[PC0 2 = 0.01] b 
AvK97[PC0 2 = 0.01] 
AVK97[PC0 2 = 0.01] 
/wK97bc0 2 = 0.01] 


0.29 
0.44 
0.80 
0.74 



a i2 is the cooling function from SMS08 and SMS09, and corresponds to present-day Earth-like greenhouse gas abundance (CO2 
partial pressure of ~3 X 10~ 4 bars). 
b AvK97 is the cooling function from Williams & Kasting (1997), and pco 2 is the partial pressure of C02- 
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Fig. 1. — Temperature evolution maps for cold-start models at 1 AU. Both models have orbital eccentricity of 0.8 along with 
Earth-like 23.5° polar obliquity and 1 bar surface pressure. Temperature is initialized to 100 K, and quickly rises to near 273 K. 
The melting of the ice-cover is handled in accordance with the prescription of §2.2. Left: CO2 partial pressure is held constant at 
0.01 bars. In this model, once the equatorial region melts, the region of surface that has melted ice-cover grows steadily until the 
entire planet has melted, and temperatures eventually grow to more than 400 K over much of the planet (not shown). Right: CO2 
partial pressure varies with temperature, in a crude simulation of a "chemical thermostat". In this model, the climate reaches a stable 
state with equatorial melt regions and polar ice-cover. 
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Fig. 2. — Compressed Milankovitch-like evolution of eccentricity and temperature at 1 AU and at 0.8 AU. Planets are initialized 
with warm equator and cold poles, similar to present-day Earth. In the top row (1 AU), the model planets are the same as in Fig. 1, 
except the eccentricity varies sinusoidally between and 0.83 with a 25-year period, to simulate a time-acceleration (by a factor of 
~10 2 to ~10 4 ) of a Milankovitch-like cycle. When the eccentricity falls below 0.05, the planet's albedo spikes to 0.8, simulating a 
catastrophic event that plunges the planet into a snowball state, with the latent heat prescription of §2.2. In the bottom row (0.8 AU), 
the eccentricity varies between 0.1 and 0.33, also with a 25-year period. Left: CO2 partial pressure is held fixed at 0.01 bars. As in 
the left panel of Fig. 1, these planets do not establish a temperate equilibrium. Right: CO2 partial pressure varies with temperature. 
Here, temperature increases are muted by reduced greenhouse effect once the ice-cover has melted somewhere. 
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Fig. 3. — Eccentricity evolution of an Earth-mass planet at 1 AU under the influence of a range of giant planet masses and orbits, 
labeled by the giant planet (GP) semimajor axes a, eccentricities e, and masses M. Top left: Effect of changing the giant planet 
mass between Saturn's mass and 3x Jupiter's mass for the case of acp = 5AU, egp = 0.4. Top right: Effect of changing the giant 
planet eccentricity between 0.05 and 0.4 for the case of oqp = 5AU, Map = Mj up . Bottom left: Two cases with similar eccentricity 
amplitudes but very different planetary system architectures, although both with Map = Mj up : acp = 0.5AU, egp = 0.1 (solid line) 
and ciqp = 5AU, ecp = 0.4 (dashed line). Bottom right: An extreme case, with acp = 30 AU.ecp = 0.925, and Map = 10 Mj up 
(dashed line), and with a = 10 AU, ecp = 0.25, Map = l0M Jup , and lap = 75° (solid line). 
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Fig. 4. — The maximum eccentricity reached by a disk of massless test particles over a million year integration of two systems 
containing a single Jupiter-mass giant planet in different configurations: acp = 5AU, egp = 0.4 (black dots), and acp = 0.4AU, egp = 
0.2 (gray dots). 



